Efficient simulations of charged colloidal dispersions: A density functional approach 
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A numerical method is presented for first-principle simulations of charged colloidal dispersions in 
electrolyte solutions. Utilizing a smoothed profile for colloid-solvent boundaries, efficient mesoscopic 
simulations are enabled for modeling dispersions of many colloidal particles exhibiting many-body 
electrostatic interactions. The validity of the method was examined for simple colloid geometries, 
and the efficiency was demonstrated by calculating stable structures of two-dimensional dispersions, 
which resulted in the formation of colloidal crystals. 
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I. INTRODUCTION 

Electrostatic interactions play a crucial role in colloidal 
dispersions Q, |^ 1^. When colloidal particles are im- 
mersed in electrolyte solutions, the so-called electric dou- 
ble layer is formed. The electric double layer is a cloud of 
counterions dissociated from the surfaces of the colloids 
into the solvent surrounding the colloidal particles. Most 
counterions are localized within the double layer due to 
the electrostatic attractive interactions between counte- 
rions and the inversely charged colloidal surfaces, while 
the entropy of the counterions tends to delocalize them. 
The thickness of the double layer is then determined by 
the competition between these conflicting effects. The 
static density profiles of the counterions can be calcu- 
lated properly using the Poisson-Boltzmann theory, and 
its linearized version leads to the well-known screened 
Coulomb interaction for a pair of likely charged colloids. 
Although the screened Coulomb potential is widely used 
to simulate charged colloidal dispersions, the lineariza- 
tion is justified only for large interparticle separations. 
Deviations from screened Coulomb interaction become 
notable for the interparticle separation smaller than the 
Debye screening length. Furthermore, many-body inter- 
actions become significant for dense colloidal dispersions. 
We thus need an alternative framework which is appli- 
cable for simulating dense dispersions composed of many 
charged colloids. 

In principle, the above problems can be resolved prop- 
erly if molecular dynamics (MD) or Monte Carlo (MC) 
simulations are used with treating counterions explicitly 
as millions of charged particles. From a computational 
point of view, however, such fully microscopic simula- 
tions are prohibitively inefficient because of the huge 
asymmetries both in size and time scales between col- 
loidal particles (large and slow) and counterions (small 
and fast). An enormous number of counterions and sim- 
ulation steps are required even for a system composed 
of only a few colloidal particles. Alternatively, counteri- 
ons can be modeled as a coarse-grained continuum object 
which is governed by a set of partial differential equations 
(PDEs) with appropriate boundary conditions defined at 
the fluid-colloid interface as demonstrated in some pre- 



vious studies ^ IE 1^ rather than microscopic particles 
governed by Newton's equations of motion. This idea can 
be most simply implemented by utilizing the finite ele- 
ment method (FEM) , which is a very natural and sensible 
method to simulate solid particles with arbitrary shapes 
in a discrete computational space. Several boundary- 
fitted unstructured mesh have been applied to specific 
problems, so that the shapes of the particles are prop- 
erly expressed I3- Thus, in principle it is possible 
to apply this method to dispersions consisting of many 
particles with any shape. However, a numerical ineffi- 
ciency arises from the following: i) re-constructions of the 
irregular mesh are necessary at every simulation step ac- 
cording to the temporal particle positions, and ii) PDEs 
must be solved under boundary conditions imposed on 
the surfaces of all colloidal particles. The computational 
demands thus are enormous for systems involving many 
particles, even if the shapes are all spherical. 

In order to overcome the problems mentioned above, 
a new idea was put forwarded by Lowcn et al. for their 
first principle MD simulations of charged colloidal disper- 
sions [iflliUl^lidjIij- Utilizing a pseudopotential, this 
method enabled us to use the conventional Cartesian co- 
ordinate and to calculate the force acting on each particle 
due to the solvents efficiently. In the present paper, we 
extend this idea by introducing a smoothed profile func- 
tion (/>(r) to represent the colloid-solvent interface with a 
finite thickness ^. The validity of the method is exam- 
ined in some simple situations with changing the interface 
thickness ^. Then the efficiency of the method is demon- 
strated by simulating two-dimensional dispersions, which 
resulted in the formation of colloidal crystals. The same 
type of method has already been applie d successfully to 
liquid-crystal colloid dispersions |15l Il6| . 

II. DENSITY FUNCTIONAL APPROACH 

A. original governing equations 

According to the density functional theory established 
for charged colloidal dispersions 0, O, E|, colloids 
are treated explicitly as particles, while counterions are 
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treated as a continuum object. Let us consider a system 
consisting N negatively charged colloidal particles with 
radius a and a solution of monovalent counterions and 
colons whose bulk concentrations are set to pq. Each col- 
loidal particles is carrying a negative charge —Ze, which 
is distributed uniformly on its surface of area A — Aira? 
(three-dimension) or — 27ra (two-dimension). Here e rep- 
resents the unit charge. The solvent is assumed to have 
an uniform dielectric constant e, and spatial distribu- 
tions of counterions(-l-) and coions(— ) are characterized 
by the local number density p+{r) and p-{r). The over- 
all charge neutrality of the system is guaranteed by the 
constraint 
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with epc{r) = e(p+(r) — p_(r)). The integral runs over 
the total volume and p±{r) = is to be strictly satisfied 
if r is inside of the particles. For a given colloidal config- 
uration {Ri, • • • , Rn}, the free energy ^ of the system 
is given by the functional of p±{r) as |3. IT^ 
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^eie = 2 I dre[pc{r) + q{r)]^{r), 
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where Tid and Teie represent the ideal gas contribution of 
the ions and the electrostatic contribution resulting from 
Coulomb interaction, respectively. The distribution of 
the surface charge of colloidal particles is given by 



Ze 



N 



{r) = -—J2Sia-\r-R,\), 



(5) 



and ip{r) is the electrostatic potential defined by the so- 
lution of the Poisson equation 



V^i^{r) = —[p^{r)+q{r)]. 



(6) 



The equilibrium density of counterions and colons 
p'±{r) are given by the solution of the variational equa- 
tion 
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where ft = !F — p J drpc{r) is the grand potential func- 
tional and p is the chemical potential. Equation ^ leads 
to 



p±ir) = po exp[(=Fei/'(r) ± p)/kBT], 



(8) 



where the constant p must be determined by the charge 
neutrality condition Eq. (Q. Substituting p±{r) for 



p±{r) in Eq. lO yields the Poisson-Boltzmann (PB) 
equation. From a computational point of view, solving 
the PB equation for dispersions of many colloidal par- 
ticles is quite demanding because the equation must be 
solved iteratively to impose correct boundary conditions 
defined on surfaces of all the colloidal particles. Usually, 
this can be done by employing non-Cartesian coordinate 
systems as in FEM 7, 8, 9], which however makes nu- 
merical calculations very complicated and inefficient for 
dispersions involving many moving particles. Another se- 
rious problem arises if one calculates the force acting on 
colloids induced by the counterions ff^ — —dfl/dRi due 
to singularities in similar to the case of liquid-crystal 
colloid dispersions [l5lll6| . 



B. smoothed profile method 

In order to improve the numerical inefficiency due to 
the moving boundary problem, a smoothed profile was in- 
troduced to the colloid-solvent interface with its thickness 
^ rather than the original sharp profile corresponds to 
^ = 0. This simple modification greatly benefits the per- 
formance of numerical computations. Similarly to earlier 
studies El El El 113, the smoothed profile function 
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is used for the individual particle i, where ^ represents 
the interface thickness. Then, the governing equations 
and lO-® can be re- written using the overall profile 
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VV(r) = -— [(1 - 0(r)) pe{r) + q{r)] . (14) 

The delta function in Eq. © is replaced with a smooth 
function |V(/)i(r)| in Eq. 1)13(1 to remove the numerical 
singularity. Since Pc{r) and q{r) are now continuous 
functions of r in the whole domain without any bound- 
ary conditions, the electrostatic potential ip{r) can be 
obtained very efficiently by using the fast Fourier trans- 
form to solve Eq. 11411 . The equilibrium density profile 
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P±{r) is calculated iteratively until Eqs. (jSJ, (|10() . and 
H14|l become self-consistent. The grand potential is given 

by 



n[p±iry,Ri,--- ,RN]=:F-fiJ dr(l 



(15) 

We note that the above equations Eqs. ((Tm) - ((Tl|l reduce 
simply to the original Eqs. 1^ and (jSJ-linii respectively, 
for C ^ 0. 

Once the equilibrium density Pj?(t") is determined, the 
force f^^ acting on the ith particle directly follows the 
Hellmann-Feynman theorem, 

dn[p''^{r);Ri,--- ,Rn] 
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One can easily calculate the force since both d(j)i{r) / dRi 
and dq{r) / dRi are analytical functions of R4. The par- 
ticle positions should be updated using the appropriate 
equations of motion. 
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where /f ^ represents the force due to direct particle- 
particle interactions, and ff and /f are the hydrody- 
namic and random forces. Repeating this procedure en- 
ables us to perform first-principle simulations for charged 
colloidal dispersions without neglecting many-body inter- 
actions. 

The present method can be regarded as an alterna- 
tive representation of the idea proposed by Lowen et al. 
\lO, 11, 12, 13, 14J, where the total charge of a single par- 
ticle is assumed to be located at the center of the particle 
and the penetration of counterions into the colloids are 
avoided effectively by using a pscudopotential between 
colloids and counterions. Kodama et al. extended this 
idea in their recent paper where the total charge 
is distributed on the colloidal surface and the pscudopo- 
tential is used to avoid the penetration. In the present 
method, the densities of counterions and colons are de- 
fined as (1 — 4i{r))p±{r) and the total charge of a colloid 
is distributed on its surface by using also the profile func- 
tion (/)(r) so that the ions cannot penetrate into colloids 
explicitly and the conservation of the ions are satisfied 
automatically. 



III. NUMERICAL CALCULATIONS 

We have performed simulations for two-dimensional 
dispersions composed of colloidal disks, which are the 




FIG. 1: Electrostatic potential tp{r) around an infinitely long 
cylindrical rod of radius a = 5 as a function of the distance 
from the center. The solid line indicates the analytic solution. 
The electrostatic potential and the length are measured in 
units of fcflT/e and As, respectively. Inset: relative error in 
the electrostatic potential for ^ = 0.5 (circles), 1.0 (squares), 
1.5 (diamonds), and 2.0 (triangles) as a function of r — a. 



two-dimensional representation of infinitely long cylin- 
drical rods. The Debye screening length is given by 
Xd — 1 / -y/STrABpo , where As = /eksT is the Bjerrum 
length and used as the unit of length in this paper. The 
radius and line charge density at the surface of the disks 
are chosen as a = 5 and Ae = Ze/Xs = e, respectively. 

To test the accuracy of the present method, we first 
consider a rod located at a center of the cylindrical 
(Wigner-Seitz) cell whose inner radius is i? = 64. Here 
we assume the system contains no additional salt. The 
analytical solution of the PB equation is known for this 
simple geometry, and the electrostatic potential is given 

by mm 
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for a < r < R. Here, the constant 7 is determined as the 
solution of the algebraic equation 
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where the constant R„i is given by 
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Furthermore, the inverse screening length kjj is defined 
by k|) = 4(1-1- 7^)/i?. The electrostatic potential is cal- 
culated using our method with ^ = 0.5, 1.0, 1.5, and 
2.0 and plotted in Fig. ^ We see that the numerical 
results are in good agreements with the analytic solu- 
tion for r — a > ^, though some deviations are found for 
r — a < ^ as an artifact of the smoothed profile. We 
emphasize that the inset of Fig. ^ shows the relative er- 
ror (-0(r,^) - V'cxact(f))/'0cxact(r) is Only within 1% for 
r — a > ^. 
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FIG. 2: The pair force f'^'^\r) between two colloidal disks 
as a function of the inter-disk separation r with ^ = 0.5 (cir- 
cles), 1.0 (squares), 1.5 (diamonds), and 2.0 (triangles). The 
dashed line indicates the LPB solution. The units of force 
and length are ksT/X^ and As, respectively. Inset: log-log 
plot of /(2)(r). 



We next consider the pair interaction between charged 
disks immersed in a solution of countcrions and colons. 
The system consists of 1024 x 1024 grid points with the 
periodic boundary condition. The linear length of the 
system is L = 512 and the Debye screening length is 
chosen as Ad = 50 in unit of A^. The equiUbrium den- 
sity profile is calculated for given inter-disk separations 
r, then, the force acting on the pair f^-^^ (r) is calculated 
using Eq. (|16|l for different values of the Interface thick- 
ness C = 0.5, 1.0, 1.5, and 2.0. In Fig. EJ the force ob- 
tained by the present method is plotted and compared 
with the analytical solution of the linearized PB (LPB) 

equation /lpb('') — —dv/dr with 
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where Ko{x) and Ki(x) are Bessel functions of imaginary 
argument. Since the thickness of the electric double layer 
is roughly given by the Debye screening length Au, the 
force obtained by our numerical method agrees well with 

(2) 

the linearized solution fi^pg{r) for r— 2a > A^i. For short 
distances r — 2a < Ad, on the other hand, deviations from 
the LPB solution become notable. It is seen in the inset 
of Fig. [3 that the deviations between numerical results 
and the LPB solution become notable also for large r. 
This is nothing more than the artifact of the periodic 
boundary condition used in our numerical calculations 
where f^'^\r) must be zero at r = L/2 — 256. 

The dependence of the interface thickness ^ on f'''^\r) 
is similar to the case of he electrostatic potential i^ir) 
shown in Fig. ^ The all numerical curves in Fig. |21 are 
almost identical (within 1% error) for r — 2a < 2f , while 
some deviations are observed for very short distances 
r — 2a < 2^. In this case, the overlapping of the in- 
terface functions 4)i{r) occurs between two disks. If we 
use a very small ^ and an infinite number of grid points. 



FIG. 3: The three-body force /'^'(/i,ri2) acting on the third 
disk as a function of h defined in the inset. Four curves are 
plotted for different values of ri2, the inter-disk separation of 
disks 1 and 2. The units of force and length are fesT/Aa and 
As, respectively. 



we may reproduce the force curve obtained by more accu- 
rate, but numerically more expensive, methods like FEM 
quantitatively for small inter-disk separations. However, 
the numerical cost would also be very expensive in such a 
case. In other words, the trade-off for the increase in nu- 
merical efficiency using the non-zero interface thickness 
is some loss of the numerical accuracy. This may give rise 
to some quantitative errors when the separation between 
colloids are very small r — 2a < 2^. 

Thirdly, we consider a system with three disks in or- 
der to examine three-body interactions acting on charged 
colloidal disks. Physical parameters were chosen identi- 
cal to the two-disk case. The geometry of the three-disk 
system is shown in the inset of Fig. O The third disk is 
located in the mid-plane of the first and the second disks 
which have fixed separations ri2 . We calculated the total 
force /3(/i, ri2) acting on the third disk with varying ft,, 
the distance from the axis connecting the first and sec- 
ond disks. Since fz{h^ri2) is always vertical because of 
the symmetry, the vertical component of the three-body 
force f'^^\h, ri2) acting on the third disk can be defined 
by 

{h, ri2) = h{h, ri2) - (ri3) cos B ~ (raa) cos 0, 

(22) 

where rij denotes the distance between the ith and jth 



disks, and 9 is defined as cos 6* ~ h/ri^. In Eq. I|22|l . 
ri3 = 7-23 and h — \/rf^ — (ryil'^Y' due to the symme- 
try of the geometry, and the pair force /^^^ is given by 
the numerical results shown in Fig. 13 Figure 13 shows 
/('^)(/i, ri2) as a function of h with four different values 
of r\2. Although we illustrate the numerical results only 
with ^ = 1.0, we have carried out the same calculations 
also with ^ = 0.5, 1.5, and 2.0. For the three-disk ge- 
ometry examined here, all the curves tend to collapse 
onto each other within 1% deviations. It is seen in Fig.|3| 
that the three-body force f^^^h,ri2) is always attrac- 
tive in this geometry. This tendency agrees well with the 
results in Ref. ^ where similar calculations have been 
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FIG. 4: Configurations of colloidal disks in the initial (a) 
and the final (b) states. The temporal charge density (1 — 
cj>{r))ep%'^{r) is shown with a color map. The scale of the 
color map is shown at the bottom in units of eAo^. 



performed by using FEM. A similar tendency has been 
observed in microscopic MD (2^J and MC ^] simula- 
tions as well as recent experiments for the same 
geometry in three-dimensional systems. 

Finally, we performed a simple demonstration for the 
crystallization of iV = 42 colloidal disks interacting each 
other. Here, the system consists of 256 x 256 grid points 
with the periodic boundary condition and the linear 
length of the system is L = 256 while the other parame- 
ters are identical to those in the previous cases of = 2 
and 3. The positions of colloidal disks are followed by 
the steepest descent-type equation of motion 



PS 



(23) 



which is obtained by substituting cPRi/dt^ — 0, /f — 
0, and — —C,dRi/dt with the friction constant C 
in Eq. 1^3). /f"^ = -dEpp/dRi is the force aris- 
ing from the potential Epp acting directly between a 
pair of colloidal disks. We defined Epp as the repul- 
sive part of the Lennard- Jones potential, Epp/kpT — 

O-^El'i' Ef=.+i[{2a/l^. - RAr- (2a/|i?. - R.lf + 
1/4] truncated at the minimum distance \Ri — Rj\ — 
2y^a. In Fig. H we show snapshots of the initial (a) 
and the final (b) configurations. Starting from a non- 
overlapping random configuration shown in Fig. E^a), 
the colloidal disks move simply to reduce the total free 
energy. Eventually, the system attains a crystalline 



state with a hexagonal close packed structure shown in 
Fig. Efb) even at a very small packing fraction rj = 
tto^N/L'^ ~ 0.05 of colloid without any effective long- 
range interactions between colloidal disks. Similar crys- 
talline structures have been observed in real experiments 
on charge-stabilized colloids . It is worth mentioning 
the computational efiiciency of our numerical method. 
The demonstration shown in Fig ^ which required 600 
time steps, takes only one hour on a PC with a single 
Pentium4 2.8GHz CPU. 



IV. CONCLUDING REMARKS 

A mesoscopic first-principle method are proposed for 
simulating charged colloidal dispersions. In order to re- 
move the numerical inefficiency due to the moving bound- 
ary condition imposed on the colloid surface, a smoothed 
profile was introduced to represent the colloid-solvent in- 
terface. In our method, the effects of counterions are 
considered within the framework of a density functional 
theory. Furthermore, many-body effects among charged 
colloids are also included properly. We have examined 
the accuracy of the present method by changing the in- 
terface thickness ^ and found that the accuracy is satis- 
factory as far as the distance between colloidal disks is 
larger than 2{a + 

Our final goal is to develop a simulation method ap- 
plicable to dynamical problems of charged colloidal dis- 
persions such as electrophoresis where the coupling be- 
tween hydrodynamics and electrostatic interactions are 
crucial In the present paper, we restricted our atten- 
tions only to static problems with employing the adia- 
batic approximation, i.e., p±{r) follows instantaneously 
to the motions of the colloidal disks or particles. This 
is OK for calculating stable colloidal structures in the 
dispersions. In the cases of dynamical problems, the 
time evolution of p±{r) should be determined by cou- 
pling equations of hydrodynamics and thermal diffusion. 
The PB equation is not appropriate for treating dynam- 
ical problems in which the counterion density becomes 
anisotropic around a particle because of the friction be- 
tween counterions and solvents. Dipoles are induced if 
this happens, and thus interactions between colloids are 
no longer screened. There appear long-range interac- 
tions between colloids, which must be important in many 
practical problems including electrophoresis for example. 
Integration of the present method and the method for 
colloids in Newtonian fluids IT^ present promising ap- 
proaches to solve these cases, and efforts to this end are 
currently underway. 

As this paper was being written for publication, we 
became aware of a parallel effort by Kodama et al. [2]| . 
While the omission of several details in their implemen- 
tation make a detailed comparison between the two ap- 
proaches impossible at this time, it will be most useful 
in the future to make a detailed comparison of the two 
methods both on efficiency and accuracy. 
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